#include<stdio.h>
#include <math.h>
#define PI 3.1415926

//返回全局变量 sunrize sunset

double sunrizetime,sunsettime;

void CalSunRise(double year,double month,double day,double jingdu,double weidu)
{
double rssj;//日升时间
double rlsj;//日落时间
double ztsj;//中天时间
double rszztsj;//日升至中天时间
double ztzrlsj;//中天至日落时间
double jsc;//均时差（小时）

double shiqu;//时区8

double JD1;//儒略日中间值
double JD2;//儒略日
double G4;//参考A
double H4;//参考B
double D8,D9,D10,D11,D12,D13,D14,D15,D16,D17,D18,D19,D20,D21,D22;
double F10,F11,F12,F13,F14,F15,F16,F17,F18,F19,F20,F21,F22;
double D25,D26,F26,H26;
double K8,K9,K10,K11,K12,K13,K14,K15,K16,K17,K18,K19,K20,K21,K22;
double M10,M11,M12,M13,M14,M15,M16,M17,M18,M19,M20,M21,M22;
double rzsjdyb;//日照时间的一半
double weiduhudu;//纬度弧度 

double xzlzxz;//修正量正弦值
double R8,R9,R10,R11,R12,R13,R14,R15,R16,R17,R18,R19,R20,R21,R22;
double T10,T11,T12,T13,T14,T15,T16,T17,T18,T19,T20,T21,T22;
//jingdu=DFM(116,23,29);//天安门
//weidu=DFM(39,54,20);
shiqu=8;

day=day -0.5;/************************************/
if(year>0)
{
year=year;
}
else
{
year=year-1;
}
if(month>2)
{
year=year;
month=month;
}
else
{
year=year-1;
month=month-12;/********************************************/
}
weiduhudu=weidu/180*PI;
G4=floor(year/100);
H4=2-G4 floor(G4/4);
JD1=floor(365.25*(year 4716)) floor(30.6001*(month 1)) day H4-1524.5;
if(JD1>=2299160.5)
{
JD2=JD1;
}
else
{
JD2=JD1-H4;
}
D8=JD2-jingdu/360;
D9=(D8-2451545)/36525;
D10=84381.448-46.815*D9-0.00059*D9*D9 0.001813*D9*D9*D9;
F10=D10/3600*PI/180;
D11=280.46645 36000.76983*D9 0.0003032*D9*D9;
F11=D11*PI/180;
D12=357.5291 35999.0503*D9 - 0.0001559*D9*D9 -0.00000048*D9*D9*D9;
F12=D12*PI/180;
D13=0.016708617 - 0.000042037*D9 - 0.0000001236*D9*D9;
F13=D13;
D14=(1.9146 - 0.004817*D9 -0.000014*D9*D9)*sin(F12) (0.019993 - 0.000101*D9)*sin(2*F12) 0.00029*sin(3*F12);
F14=D14*PI/180;
D15=D11 D14;
F15=D15*PI/180;
D16=D12 D14;
F16=D16*PI/180;
D17= 1.000001018*(1-F13*F13) / (1 F13*cos(F16));
F17=D17;
D18= 125.04 - 1934.136*D9;
F18=D18*PI/180;
D19=D15 - 0.00569 -0.00478*sin(D18);
F19=D19*PI/180;
F20=asin(sin(F10)*sin(F15));
D20=F20*180/PI;
D21=D10/3600 0.00256*cos(F18);
F21=D21/180*PI;
F22=asin(sin(F21)*sin(F19));
D22=F21/PI*180;
D25=(tan(F10/2))*(tan(F10/2));//sin/cos
D26=D25*sin(2*F11)-2*F13*sin(F12) 4*F13*D25*sin(F12)*cos(2*F11)-(D25*D25/2)*sin(4*F11)-(5/4*F13*F13)*sin(2*F12);
F26=D26*180/PI;
H26=F26*4;
//================
//xzlzxz=sin((16/F17 34)*-1/60*PI/180)；//34大气折射修正 单位：分
xzlzxz=sin((16/F17 34)*-1/60*PI/180);
rzsjdyb=acos((xzlzxz-sin(F22)*sin(weiduhudu))/(cos(F22)*cos(weiduhudu)))/PI*12;
K8=D8-rzsjdyb/24;
K9=(K8-2451545)/36525;
K10=84381.448-46.815*K9-0.00059*K9*K9 0.001813*K9*K9*K9;
M10=K10/3600*PI/180;
K11=280.46645 36000.76983*K9 0.0003032*K9*K9;
M11=K11*PI/180;
K12=357.5291 35999.0503*K9 - 0.0001559*K9*K9 -0.00000048*K9*K9*K9;
M12=K12*PI/180;
K13=0.016708617 - 0.000042037*K9 - 0.0000001236*K9*K9;
M13=K13;
K14=(1.9146 - 0.004817*K9 -0.000014*K9*K9)*sin(M12) (0.019993 - 0.000101*K9)*sin(2*M12) 0.00029*sin(3*M12);
M14=K14*PI/180;
K15=K11 K14;
M15=K15*PI/180;
K16=K12 K14;
M16=K16*PI/180;
K17= 1.000001018*(1-M13*M13) / (1 M13*cos(M16));
M17=K17;
K18= 125.04 - 1934.136*K9;
M18=K18*PI/180;
K19=K15 - 0.00569 -0.00478*sin(K18);
M19=K19*PI/180;
M20==asin(sin(M10)*sin(M15));
K20=M20*180/PI;
K21=K10/3600 0.00256*cos(M18);
M21=K21/180*PI;
M22=asin(sin(M21)*sin(M19));
K22=M22/PI*180;
//============================================
R8=D8 rzsjdyb/24;
R9=(R8-2451545)/36525;
R10=84381.448-46.815*R9-0.00059*R9*R9 0.001813*R9*R9*R9;
T10=R10/3600*PI/180;
R11=280.46645 36000.76983*R9 0.0003032*R9*R9;
T11=R11*PI/180;
R12=357.5291 35999.0503*R9 - 0.0001559*R9*R9 -0.00000048*R9*R9*R9;
T12=R12*PI/180;
R13=0.016708617 - 0.000042037*R9 - 0.0000001236*R9*R9;
T13=R13;
R14=(1.9146 - 0.004817*R9 -0.000014*R9*R9)*sin(T12) (0.019993 - 0.000101*R9)*sin(2*T12) 0.00029*sin(3*T12);
T14=R14*PI/180;
R15=R11 R14;
T15=R15*PI/180;
R16=R12 R14;
T16=R16*PI/180;
R17= 1.000001018*(1-T13*T13) / (1 T13*cos(T16));
T17=R17;
R18= 125.04 - 1934.136*R9;
T18=R18*PI/180;
R19=R15 - 0.00569 -0.00478*sin(R18);
T19=R19*PI/180;
T20=asin(sin(T10)*sin(T15));
R20=T20/PI*180;
R21=R10/3600 0.00256*cos(T18);
T21=R21/180*PI;
T22=asin(sin(T21)*sin(T19));
R22=T22/180*PI;
rszztsj=acos((xzlzxz-sin(M22)*sin(weiduhudu))/(cos(M22)*cos(weiduhudu)))/PI*12;
ztzrlsj=acos((xzlzxz-sin(T22)*sin(weiduhudu))/(cos(T22)*cos(weiduhudu)))/PI*12;
//? jsc;//均时差（小时）
jsc=H26/60;
//? ztsj;//中天时间
ztsj=12-jsc-(jingdu-shiqu*15)*4/60;
rssj=ztsj-rszztsj;//日升时间=中天时间-日升至中天时间
rlsj=ztsj ztzrlsj;//日落时间=中天时间 中天至日落时间
sunrizetime=rssj;
sunsettime=rlsj;
}

void printtime(double t)
{
printf("%d:%d:%d\n",(int)(floor(t)),(int)floor((t-floor(t))*60),(int)floor((t*60-floor(t*60))*60));
}

double DFM(double d,double f,double m)
{
return d f/60 m/3600;
}



void main()
{
CalSunRise(2013,11,11,DFM(116,23,29),DFM(39,54,20));
printtime(sunrizetime);
printtime(sunsettime);
} 
